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Abstract 

Locally refined meshes impose severe stability constraints on explicit time-stepping methods for the numer- 
ical simulation of time dependent wave phenomena. Local time-stepping methods overcome that bottleneck 
by using smaller time-steps precisely where the smallest elements in the mesh are located. Starting from 
classical Adams-Bashforth multi-step methods, local time-stepping methods of arbitrarily high order of ac- 
curacy are derived for damped wave equations. When combined with a finite clement discretization in 
space with an essentially diagonal mass matrix, the resulting time-marching schemes are fully explicit and 
thus inherently parallel. Numerical experiments with continuous and discontinuous Galerkin finite element 
discretizations validate the theory and illustrate the usefulness of these local time-stepping methods. 
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1. Introduction 

Efficient numerical methods for the simulation of damped wave phenomena are of fundamental impor- 
tance in acoustic, electromagnetic or seismic wave propagation. In the presence of complex geometry, such 
as cracks, sharp corners or irregular material interfaces, standard finite difference methods generally become 
ineffective and cumbersome. In contrast, finite element methods (FEMs) easily handle unstructured meshes 
and local refinement. Moreover, their extension to high order is straightforward, a key feature to keep 
numerical dispersion minimal. 

The finite element Galerkin approximation of hyperbolic problems typically leads to a system of ordinary 
differential equations. However, if explicit time-stepping is subsequently employed, the mass matrix arising 
from the spatial discretization by standard conforming finite elements must be inverted at each time-step: a 
major drawback in terms of efficiency. To overcome that difficulty, various "mass lumping" techniques have 
been proposed, which effectively replace the mass matrix by a diagonal approximation. While straightfor- 
ward for piecewise linear elements [TJ [5] , mass lumping techniques require special quadrature rules at higher 
order to preserve the accuracy and guarantee numerical stability 0] . 

Discontinuous Galerkin (DG) methods offer an attractive and increasingly popular alternative for the 
spatial discretization of time-dependent hyperbolic problems [51 [51 [71 [HI [HI [TO]- Not only do they accommodate 
elements of various types and shapes, irregular non-matching grids, and even locally varying polynomial 
order, and hence offer greater flexibility in the mesh design. They also lead to a block-diagonal mass 
matrix, with block size equal to the number of degrees of freedom per element; in fact, for a judicious 
choice of (locally orthogonal) shape functions, the mass matrix is truly diagonal. Thus, when a spatial DG 
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discretization is combined with explicit time integration, the resulting time-marching scheme will be truly 
explicit and inherently parallel. 

In the presence of complex geometry, adaptivity and mesh refinement are certainly key for the efficient 
numerical simulation of wave phenomena. However, locally refined meshes impose severe stability constraints 
on explicit time-marching schemes, where the maximal time-step allowed by the CFL condition is dictated 
by the smallest elements in the mesh. When mesh refinement is restricted to a small region, the use of 
implicit methods, or a very small time-step in the entire computational domain, are very high a price to pay. 
To overcome this overly restricitve stability constraint, various local time-stepping (LTS) schemes pTT l [T2 | IT3 ) 
were developed, which use either implicit time-stepping or explicit smaller time-steps, but only where the 
smallest elements in the mesh are located. 

Since DG methods are inherently local, they are particularly well-suited for the development of explicit 
local time-stepping schemes [TO]. By combining the sympletic Stormer-Verlet method with a DG discretiza- 
tion, Piperno derived a symplectic LTS scheme for Maxwell's equations in a non-conducting medium [14], 
which is explicit and second-order accurate. In [15J . Montseny et al. combined a similar recursive inte- 
grator with discontinuous hexahedral elements. Starting from the so-called arbitrary high-order derivatives 
(ADER) DG approach, alternative explicit LTS methods for Maxwell's equations [TB] and for elastic wave 
equations [T7] were proposed. In [TS] , the LTS approach from Collino et al. [TT] [T^] was combined with a 
DG-FE discretization for the numerical solution of symmetric first-order hyperbolic systems. Based on en- 
ergy conservation, that LTS approach is second-order and explicit inside the coarse and the fine mesh; at the 
interface, however, it nonetheless requires at every time-step the solution of a linear system. More recently, 
Constantinescu and Sandu devised multirate explicit methods for hyperbolic conservation laws, which are 
based on both Runge-Kutta and Adams-Bashforth schemes combined with a finite volume discretization 
[Hfl |2"U] . Again these multirate schemes are limited to second-order accuracy. 

Starting from the standard leap-frog method, Diaz and Grote proposed energy conserving fully explicit 
LTS integrators of arbitrarily high accuracy for the classical wave equation [5T] ; that approach was extended 
to Maxwell's equations in [32] for non- conductive media. By blending the leap-frog and the Crank-Nicolson 
methods, a second-order LTS scheme was also derived there for (damped) electromagnetic waves in con- 
ducting media, yet this approach cannot be readily extended beyond order two. 

To achieve arbitrarily high accuracy in the presence of dissipation, while remaining fully explicit, we shall 
derive here explicit LTS methods for damped wave equations based on Adams-Bashforth (AB) multi-step 
schemes. The rest of the paper is organized as follows. In section 2, we first recall the standard continuous 
and the symmetric interior penalty (IP) DG finite element discretizations of the second-order damped wave 
equation; there, we also rewrite the damped wave equation as a first-order hyperbolic system and recall its 
nodal DG formulation. Starting from the Adams-Bashforth multi-step schemes, we then derive in Section 3 
a LTS approach of arbitrarily high accuracy. When combined with a finite element discretization in space 
with a diagonal mass matrix, the resulting time-marching schemes remain fully explicit. Finally in Section 4, 
we present numerical experiments in one and two space dimensions, which validate the theory and underpin 
both the stability properties and the usefulness of these Adams-Bashforth LTS schemes. 

2. Finite element discretizations for the damped wave equation 

We consider the scalar damped wave equation 

utt + au t - V • (c 2 Vu) = / in Ox (0, T) , 

u(;t) = on 3!!x (0,T) , (1) 
w(-,0) = u , u t (-,0) = v infi, 

where SI is a bounded domain in M. d , d = 1,2,3. Here, / £ L 2 (0,T; L 2 (il)) is a (known) source term, 
while wo £ i?o(Q) and vq £ L 2 (£l) are prescribed initial conditions. At the boundary, d£l, we impose a 
homogeneous Dirichlct boundary condition, for simplicity. The damping coefficient, a = a(x), is assumed 
non- negative (a > 0) whereas the speed of propagation, c = c(x), is piecewise smooth and strictly positive 
(c(ar) > c > 0). 
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We shall now discretize ([I]) in space by using any one of the following three distinct FE discretizations: 
continuous (i^-conforming) hnite elements with mass lumping, a symmetric IP-DG discretization, or a 
nodal DG method. Thus, we consider shape-regular meshes Th that partition the domain fl into disjoint 
elements K, such that = UxeTh-K- The elements are triangles or quadrilaterals in two space dimensions, 
and tetrahedra or hexahedra in three dimensions, respectively. The diameter of element K is denoted by 
1%K and the mesh size, h, is given by h = max^GT^ hfc- 

2.1. Continuous Galerkin formulation 

The continuous (H ^conforming) Galerkin formulation of ([T| starts from its weak formulation: find 
u e [0,T] H£(tt) such that 

(u tt) <p) + (*u t) <p) + (cVu,cV<p) = (f,<p) V<p€H%(Sl), t e (0,T) , 
u(-,0) = u , u t (-,0) = w , 

where (•, •) denotes the standard L 2 inner product over £1. It is well-known that Q is well-posed and has a 
unique solution |23j . 

For a given partition Th of fi, assumed polygonal for simplicity, and an approximation order t > 1, we 
shall approximate the solution u(-,t) of |2| in the finite element space 

V h := {<p € Hfcfl) : <p\ K eS e (K) V K e %} , 

where S e (K) ist the space V l {K) (for triangles or tetrahedra) or Q l (K) (for quadrilaterals or hexahedra). 
Here, we consider the following semi-discrete Galerkin approximation of (pi): find u h : [0,T] — > U h such that 

(^) + (™^) + (cW\cV^) = (/,^) Vy>ev\ te(o,T), 
u h (-,o) = n /l «o ) u£(-,o) = iw 

Here, denotes the L 2 -projection onto V h . 

The semi-discrete formulation (|3| is equivalent to the second-order system of ordinary differential equa- 
tions 

M~(i)+M a ~(t) + KU(t) = F(t), te(0,T), 

MU(0) = «S, M ^ I(0) = U «- 

Here, U denotes the vector whose components are the coefficients of u h with respect to the finite element 
basis of Vh, M the mass matrix, K the stiffness matrix, whereas M„ denotes the mass matrix with weight 
(7. The matrix M is sparse, symmetric and positive definite, whereas the matrices K and M CT are sparse, 
symmetric and, in general, only positive semi-definite. 

Usually, the mass matrix M is not diagonal, yet needs to be inverted at every time-step of any explicit 
time integration scheme. To overcome this diffculty, various mass lumping techniques have been developed 
[2l[24j[25j[26], which essentially replace M with a diagonal approximation by computing the required integrals 
over each element K with judicious quadrature rules that do not effect the spatial accuracy |27) . 

2.2. Interior penalty discontinuous Galerkin formulation 

Following [S] we briefly recall the symmetric interior penalty (IP) DG formulation of (JlJ. For simplicity, 
we assume in this section that the elements are triangles or parallelograms in two space dimensions and 
tetrahedra or parallelepipeds in three dimensions, respectively. Generally, we allow for irregular (fc-irregular) 
meshes with hanging nodes [5H]. We denote by £^ the set of all interior edges of Th, by £® the set of all 
boundary edges of Th, and set Eh = U . Here, we generically refer to any element of £h as an "edge", 
both in two and three space dimensions. 
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For a piecewise smooth function ip, we introduce the following trace operators. Let e £ £^ be an interior 
edge shared by two elements K + and K~ with unit outward normal vectors , respectively. Denoting by 
v the trace of v on dK ± taken from within , we define the jump and the average on e by 

[<p] := <p+n+ + <p~n- , {<p} := {<p+ + <p~)/2 . 

On every boundary edge e £ £%, we set {ipj := <^n and ^[<p~fy := Here, n is the outward unit normal to 
the domain boundary 9f2. 

For a piecewise smooth vector- valued function 0, we analogously define the average across interior faces 
by fV'J' := + an d on boundary faces we set f := V'- The jump of a vector-valued function 

will not be used. For a vector- valued function with continuous normal components across a face e £ Eh-, 
the trace identity 

tp + (n+ • -0+) + ip~ (n- • -0") - M ■ {{0| on e , 

immediately follows from the above definitions. 

For a given partition Th of O and an approximation order £ > 1, we wish to approximate the solution 
u(t, •) of in the finite element space 

V h := {^ G L 2 («) : p|jr G S e (K) VK £ %} , 

where S (K) is the space V (K) of polynomials of total degree at most i on K if if is a triangle or 
a tetrahedra, or the space Q e (K) of polynomials of degree at most I in each variable on if if if is a 
parallelogram or a parallelepiped. Thus, we consider the following (semidiscrete) DG approximation of Q: 
find u h : [0, T] -> V' 1 such that 

(t&¥0 + (at^)+a A (u fc , ?) = (/,¥>) V^GV h , i€(0,T), 

( 4 ) 

u ft (-,o) = n fc «o, u?(-,o) = rw 

Here, again denotes the L 2 -projection onto V h whereas the DG bilinear form a^(-, •), defined on V h x V h , 
is given by 



a h (u,cp) := / c 2 Vu • V95 rfx - ^ / [u] • {{c 2 V<^J 



dA 

(5) 



e££ h e ee£ h 

The last three terms in |5]) correspond to jump and flux terms at element boundaries; they vanish when 
u, ip, £ Hq(Q) n H 1+m (Q) for m > |. Hence, the above semi-discrete DG formulation Q is consistent with 
the original continuous problem ([2]). 

In ([5| the function a penalizes the jumps of u and v over the faces of Th- To define it, we first introduce 
the functions h and c by 

( mm{h K +,h K -}, e £ £f , ( max{c\ K +(x),c\ K -(x)}, e £ £% , 

{ h K , e££^, { c\ K {x), e££f l . 

Then, on each e G £h, we set 

a| e — ac 2 !!- 1 , (6) 

where a is a positive parameter independent of the local mesh sizes and the coefficient c. There exists a 
threshold value a m i n > 0, which depends only on the shape regularity of the mesh and the approximation 
order t such that for a > a m i n the DG bilinear form ah is coercive and, hence, the discretization stable 
p9l 130] . Throughout the rest of the paper we shall assume that a > a m i n s ° that the semi-discrete problem 
Q has a unique solution which converges with optimal order [5J |SJ [3TJ 132] • In [HI 132] , a detailed convergence 



analysis and numerical study of the IP-DG method for ([5| with a — was presented. In particular, optimal 
a-priori estimates in a DG-energy norm and the L 2 -norm were derived. This theory immediately generalizes 
to the case a > 0. For sufficiently smooth solutions, the IP-DG method |5j thus yields the optimal L 2 -error 
estimate of order <D(h l+1 ). 

The semi-discrete IP-DG formulation Q is equivalent to the second-order system of ordinary differential 
equations 

M^(f)|M^(i)+KU(() = F(i), i€ (0,T), 
MU(0)=4, M^(0) = <. 

Again, the mass matrix M is sparse, symmetric and positive definite. Yet because individual elements 
decouple, M (and M^) is block-diagonal with block size equal to the number of degrees of freedom per 
element. Thus, M can be inverted at very low computational cost. In fact, for a judicious choice of (locally 
orthogonal) shape functions, M is truly diagonal. 

2.3. Nodal discontinuous Galerkin formulation 

Finally, we briefly recall the nodal discontinuous Galerkin formulation from |10j for the spatial discretiza- 
tion of ([!]) rewritten as a first-order system. To do so, we first let v :— u t , w := — V«, and thus we rewrite 
(JlJ as the first-order hyperbolic system: 

v t + av + V • (c 2 w) = / in ft x (0, T) , 

w t + V« = in fix (0,T), 

v{-,t) = 0, w(-,t) = on<9fix(0,T), ^ ' 

v(-,0) = v , w(-,Q) = -Vu infi, 

or in more compact notation as 

qi + Eq + V-J(q) = S, (8) 

with q = (u,w) T . Following [10] . we now consider the following nodal DG formulation of find q h : 
[0, T] -> V' 1 such that 

(qJ l ,^) + (Sq h ,^)+^(q h ,t/;) = (S,t/>) Vt/> G V\ t g (0, T) . (9) 

Here ~V h denotes the finite element space 

V' 1 := {V> e L 2 (n) d+1 : il>\ K G S e (K) d+1 VK e %} 

for a given partition 7ft, of f2 and an approximation order I > 1. The nodal-DG bilinear form aft(-, ■) is 
defined on V' 1 x as 

5ft(q, </>):= I ( V - H^))"4>dx- / (n ■ J-(q) - (n • ^(q))*)-^^, 

K£Th JK ee£ h Je 

where (n • .F(q))* is a suitably chosen numerical flux in the unit normal direction n. The semi-discrete 
problem |9| has a unique solution, which converges with optimal order in the L 2 -norm [10] . 

The semi-discrete nodal DG formulation ^ is equivalent to the first-order system of ordinary differential 
equations 

M^(() + M,Q(i) + CQ(f) = F(f), t€(0,T). (10) 

Here Q denotes the vector whose components are the coefficients of q h with respect to the finite element 
basis of ~V h and C the DG stiffness matrix. Because the individual elements decouple, the mass matrices M 
and M CT are sparse, symmetric, positive semi-definite and block-diagonal; moreover, M is positive definite 
and can be inverted at very low computational cost. 
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3. High-order explicit local time-stepping 



The H 1 -conforming and the IP-DG finite element discretizations of Q presented in Section 2 lead to 
the second-order system of differential equations 



M # (t) + M '¥ (i) + KU(i)=(l ' *e(o,T), 

whereas the nodal DG discretization leads to the first-order system of differential equations 



M=£(t) + M CT Q(i) + C Q(i) = . 
at 



te (0,T) 



(11) 



(12) 



In both (11) and (12) the mass matrix M is symmetric, positive definite and essentially diagonal; thus. 



M _1 or MT 2 can be computed explicitly at a negligible cost; for simplicity, we restrict ourselves here to 
the homogeneous case, i.e. F(t) = 0. 

If we multiply ( 11 1 by M~2 , we obtain 



^(t)+D-(t) + A.(t) = i 



(13) 



with z(i) = MiU(t), D = M~2M CT M~2 and A = M~3KM~2. Note that A is also sparse and symmetric 
positive semi-definite. Thus, we can rewrite (131 as a first-order problem of the form 



dy 

dt 



with 



y(t)=[z(t),-(t) 



(*) = By(i) 



B 



(14) 



I 
A D 



Similarly, we can also rewrite (12) as in the form (14) with y(<) = Q(i) and B = M 1 (— M CT — C). Hence 
all three distinct finite element discretizations from Section 2 lead to a semi-discrete system as in ( 14 ) . 



Starting from explicit multi-step Adams-Bashforth methods, we shall now derive explicit local time- 
stepping schemes of arbitrarily high accuracy for a general problem of the form 



dy 

dt 



(t) = By(t) : 



te (0,T). 



(15) 



3.1. Adams-Bashforth methods 

First, we briefly recall the construction of the classical fc-step (fcth-order) Adams-Bashforth method for 
the numerical solution of (15) [33]. Let ti = iAt and y„, y„_i,..., y ?1 -fc+i the numerical approximations to 



the exact solution y{t n ),..., y(< n -fc+i)- The solution of (15) satisfies 



y(t n + £At) = y(t n ) + ^ By(«) dt , 



<£< 1. 



(16) 



We now replace the unknown solution y(t) under the integral in ( 16 ) by the interpolation polynomial p(t) 
through the points (ti, y,), i = n — k + 1, . . . , n. It is explicitly given in terms of backward differences 



VV = y„ , V J+1 yn = V J y n - V*y»-i 



by 



k-l 



3=0 

6 



n 3 \ V3y n 



Integration of ( 16 1 with y(t) replaced by p(t) then yields the approximation y n +£ of y(t n +£At), < £ < 1, 



k-l 



y«+c = y» + AfB ^ 7i(e)V J y» , 

3=0 



(17) 



where the polynomials 7j(^) are defined as 



7i(0 = 



They are given in Table [T] for j < 3. After expressing the backward differences in terms of y„_j and setting 
£ = 1 in (17), we recover the common form of the fc-step Adams-Bashforth scheme [33J 



fc-i 



y n +i = y n + AtB 2J aj-y n - 



(18) 



where the coefficients otj, j — 0, ...,k — 1 for the second, third- and fourth-order (fc = 2,3,4) Adams- 
Bashforth schemes are given in Table [5] For higher values of k we refer to [33] . 
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Table 1: Coefficients 7j(£) for the explicit Adams-Bashforth methods. 





a 


«1 


a t 


a 3 


k 


= 2 


3 
2 


1 

2 








k 


= 3 


23 
12 


16 
12 


5 

12 





k 


= 4 


55 
24 


59 
24 


37 
24 


9 

24 



Table 2: Coefficients for the fc-th order Adams-Bashforth methods. 



3.2. Adams-Bashforth based LTS 

Starting from the classical Adams-Bashforth methods from Section 3.1, we shall now derive LTS schemes 



of arbitrarily high accuracy for (15), which allow arbitrarily small time-steps precisely where small elements 



in the spatial mesh are located. To do so, we first split the unknown vector y(t) in two parts 

y(t) = (I - P)y(t) + Py(f) = y [coarscl (f) + y [finc] (t) , 

where the matrix P is diagonal. Its diagonal entries, equal to zero or one, identify the unknowns associated 
with the locally refined region, where smaller time-steps are needed. Hence P corresponds to a discrete 
partition of unity of the degrees of freedom associated with Vh ■ 



The exact solution of ( 15 ) again satisfies 



/■In -1-5*11 

y(t n + (At) = y(t n ) + Jf B (y[<™] (t) + y N ( t ) j dt , < g < 1 . (19) 
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3 


12 3 


7i(0 





Table 3: The polynomial coefficients 7j(§) 



Since we wish to use the standard fc-step Adams-Bashforth method in the coarse region, we approximate 
the term in (191 that involve y[ coarsc l(t) as in (161, which yields 



fc-i 



y(t„ + £Ai) « y„ + A*B(I - P) 5>;(£)V J y r , 

4=0 



t„+£At 



BPy(f) di . 



(20) 



To circumvent the severe stability constraint due to the smallest elements associated with y^ finc '(t), we shall 
now treat y[ fine l(t) differently from y[ coarse l(t) and instead we approximate the integrand in (20) as 



t„+£At 



BPy(i) dt 



{At 



BPy(r)dr, 



where y(r) solves the differential equation 

dy 



k-l 



(r) = B(I - P)^7J (^) VVn + BPy(r) , 



dr 

y(0) = y„ , 



with coefficients 



^(0 = ^(0 = 1 (-if 



o v J 



; rf S =(-iy 



The polynomials 7j(£) are given in Table [| for j < 3. Replacing y(t) by y(t) in ([^Ojl, we obtain 

fe-l /-{At 
y(t„ + £At) w y„ + At B(I - P) V 7j (£) V J y„ + / BPy(r) dr . 



By considering (21 ) in integrated form, we find that 

fc-l / r £A* 



y(£At)=y(0) + B(I-P)^ U ^ (— J dr j V'y n + ^ BPy(r)rfr 



fc-i 



k-l /.{At 

y n + AtB(I-P)V 7i (£)V J y n + / BPy(r)dr. 



^From the comparison of ( 23 1 and ( 24 1 we infer that 

y(t„ + £At) « y(£At) 



(21) 



(22) 



(23) 



(24) 



Thus to advance y(t n ) from t„ to i„ + At, we shall evaluate y(At) by solving (21 ) on [0, At] numerically. We 
solve (21 ) until t — At again with a /c-step Adams-Bashforth scheme, using a smaller time-step Ar = At/p, 
where p denotes the ratio of local refinement. For m = 0, . . . ,p — 1 we then have 



k-l k-l 

y(m+i)/p = ym/p + At B(I - P) ^2 a t ^2 li 



m — i 
P 



k-l 



V-^ + ArBP^T 



l)/p: 



(25) 



where ct£, I = 0, . . . , k — 1 denote the coefficients of the classical fc-step Adams-Bashforth scheme (see Table 
p2|). Finally, after expressing the backward differences in terms of y n -i, we find 



fe-i 



fe-i 



y(m+l)/p Jm/p i 

At B(I - P) Pm,e y n -t + Ar BP ^ a e y (m _ ;)/p , 
1=0 e=o 

where the constant coefficients f3 m ,l, m = 0, . . . ,p — 1, I = 0, . . . , k — 1, satisfy 

fc-i fc-i 



m = 0, 



,P-1, 



(26) 



i=o j=e 



i I J 

£ 



1i 



m 



P 



(27) 



with 7j defined in (22). 

In summary, the LTS-ABfc(p) algorithm computes y n+ i ~ y(t„ + At), given y„, y n _ lr .., y„_& +1 , 
B (! - P)y«-i B(I - P)y„-fe+i and Py„_i/ P , Py n _ 2 / P ,---, Pyn-(fc-i)/ P as follows: 

LTS-ABfc(p) Algorithm 

1. Set y := Yn, Y-e/ P ■= ^Yn-i/p, t = 1, ■ ■ ■ ,k - 1. 

2. Setw n - e := B(l - P)y n _i, £ = 

3. Compute w„ := B(I — P)y n - 

4. For m = 0, . . . ,p — 1, compute 



l,...,k-l. 



y(m+l)/p 



Af 



fe-1 



ym/p H > , Pm.l W n _ 

^ £=0 



Af 



fe-1 



BP^a £ y ( 



m—l)/p • 



1=0 



5. Sef y„ + i := y : . 



Steps 1-4 correspond to the numerical solution of (21 1 until t — At with the fc-step Adams-Bashforth 
scheme, using the local time-step Ar = At/p. For P = or p = 1, that is without any local time-stepping, 
we thus recover the standard k-step Adams-Bashforth scheme. If the fraction of nonzero entries in P is small, 
the overall cost is dominated by the computation of w n in Step 3, which requires one multiplications by 
B(I — P) per time-step At. All further matrix- vector multiplications by BP only affect those unknowns that 
lie inside the refined region, or immediately next to it; hence, their computational cost remains negligible 
as long as the locally refined region contains a small part of ft. 

3.3. Examples of LTS-ABk(p) schemes 

We have shown above how to derive LTS-ABfc(p) schemes of arbitrarily high accuracy. Since the third- 
and fourth-order LTS-ABfc(p) schemes are probably the most relevant for applications, we now describe in 
detail the LTS-ABfc(p) schemes for k — 3, 4 and p = 2. 

For k = 3 and p = 2, we find from (251 and (26) for the first half-step of size Ar = At/2 that 

yi/2 = yo + ^B(I - P) [ao (7o(0)V°y» + 7i(0)V 1 yn + 7 2 (0)V 2 y„) 

+ ai (7o(-l/2)V°y« + ^x{-l/2)V x y n + 7 2 (-l/2)V 2 y„) 
+ « 2 (7o(-l)V°y n + 7i(-l)VVn + 7 2 (-l)V 2 y„ 



At 

At 



BP 



«oyo + aiy-1/2 + a> 2 y-i 



B(i - P) /3 , y» + /3o,iy„-i + /3o, 2 yn- 



At 



BP 



aoy« + aiy-1/2 + a 2 yn-i 
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tn+1 



re+l/2 


Yn+1 


yi 








71 




yi/2 








n-1/2 


y n 


yo 








t»-l 




y-1/2 








'n-2 


y„-i 
y„.-2 


y-i 









X 



Figure 1: Illustration of the LTS-AB3(2) scheme in one space dimension. 



with 



A),o - "o {70 (0) +7i (0) +72(0)} + ai {^(-1/2) +71 (-1/2) + ^(-1/2)} 
+ a 2 {7o(-l) +7i(-l) +72(-l)} 

1 1 



l( 1 + o + o,-H 



1 



8 



A>,1 = «o {-71 (0) - 272(0)} + ai {-7i(-l/2) - 27 2 (-l/2)} + a 2 {-%(-!) - ^(-1)} 



^{o + o}- 16 !^ 1 

12 1 ' 12 1 2 4 



s< 1+0 >— ^ 



~ , s ~ , , s ~ , * 23 16 f 11 5 2 
A>,2 = «o7 2 (0) + ai 72 (-l/2) + a 272 (-l) = ^- °- 1 ^\si + 12" 0= 12" 

Next, we perform a second half-step thus completing the full step of size At: 
_ At 



yi = yi/2 



B(I - P) |o<, (7o(l/2)V°y« + 7i(l/2)V 1 yn + 72(l/2)V 2 y„) 
+ 01 (70(0) V°yn + 7i(0)V 1 y, l +72(0)V 2 y«) 
+ a 2 (7o(-l/2)V°y» + 7i(-l/2)V 1 y« + 7 2 (-l/2)V 2 y„) 



yi/2 



— BP 

2 

At 



«oyi/2 + "iyo + Q!2y-x/2 



B(I-P) 



PlflYn + /3i,iy n -i + PlfiYn-l 



At. 



-BP 



aoyi/2 + «iyn + a 2 y-i/2 



with 



ft.o = "o {7o(l/2) + 7! (1/2) + 72(1/2)} + a x {70(0) + 71 (0) + 7 2 (0)} 
+ a 2 {7o(-l/2) + 71 (-1/2) + 72(-l/2)} 



23 
12 



1 1 

'-2-8 



29 
12 ' 



a {-7i(l/2) - 272(1/2)} + a x {-71(0) - 2^(0)} + a 2 {-^(-1/2) - 2j 2 (-l/2)} 



23 
12 



1 3 

2 ~ 4 



16 r ,5 

12<°-°> + 12 



^1,2 = a 7 2 (l/2) + Qi7 2 (0) + a 2 7 2 (-l/2) 



1 1 

2 + 4 
23 

~ 12 
10 



25 
12' 
16 
' 12 ' 



5 

12 



12 



Recall that the coefficients at correspond to the standard coefficients of the Adams-Bashforth methods given 
in Table [2j After simplification, the LTS-AB3(2) method then reads (see Figure [lj: 



yi/2 



Yn+i — yi — yi/2 



At 



-B(I-P) 



17 



12 



Yn-l 



12 



y«-2 



At 



BP 



23 , 
V2 : 



16. 
12 



At 



B(I-P) 



29 



25 



22^™ 12^™ 1 i2"^ n— 2 



At. 



BP 



23. 
12 



yi/2 



y-1/2 - 

16 

12 Yn 



For the case with k = 4 and p — 2, similar calculations yield the LTS-AB4(2) scheme: 



yi/2 = y-n 



yn+i — yi — yi/2 



At 

~2~ 

At 

~2~ 

At 

~2 

At 



B(I-P) 

55 



297 



BP 



24 

B(I-P) 



192 
59 
24 



187 
192 



Yn-l 



107 
192 



Yri-2 



25 
192 



y«-3 



y-1/2 



37 

24 



y n -i 



24 



y-3/2 



583 
192 



y n 



757 
192 



Yn-l 



485 
192 



y n -2 



119 
192 



y™-3 



BP 



24 



yi/2 



59 

24 



37, 



Yn + ^y-1/2 



24 



Yn-l 



12 Yn-l 
5 _ 

^y-1/2 



Other examples of LTS Adams-Bashforth schemes are listed in Table |4j where the coefficients of the 
schemes are given for different values of k and p. 



Scheme 


m 


Pm,0 


/3m, 1 


Pm,2 


ftm, 3 


LTS-AB2(2) 





5 

4 


1 

4 










1 


7 

4 


3 

4 








LTS-AB2(3) 





7 
6 


1 

6 










1 


9 
6 


3 
6 










2 


11 

6 


5 
6 








LTS-AB3(2) 





17 
12 


7 
12 


2 
12 







1 


29 
12 


25 
12 


8 

12 





LTS-AB3(3) 





137 
108 


40 
108 


11 
108 







1 


203 
108 


136 
108 


41 
108 







2 


281 
108 


256 
108 


83 
108 





LTS-AB4(2) 





297 
192 


187 
192 


107 
192 


25 
192 




1 


583 
192 


757 
192 


485 
192 


119 
192 


LTS-AB4(3) 





871 
648 


387 
648 


213 
648 


49 
648 




1 


1425 
648 


1437 
648 


867 
648 


207 

648 




2 


2159 
648 


2955 
648 


1917 

648 


473 

648 



Table 4: Coefficients of LTS-ABfc(p) schemes for k = 2,3,4 and p = 2, 3, 4. 



11 



3.4- Accuracy of the LTS scheme 

We now establish the accuracy of the LTS-ABfc(p) scheme. To do so, we first recall that for k > 1, we 
have 

fc-i 

J2 a i = 1 ' 



(28) 



where etj, j = 0, . . . , k — 1 are the standard coefficients of the fc-step Adams-Bashforth scheme (18); see 
Theorem 2.4 in [33J for details. Next, we shall need the following identity. 



Lemma 1. For k > 1 and p > 2 we have 



p-i 



^ P m ,£ =pae, i = 0,...,k-l, 



(29) 



m— 



where ct£ and fi m j, (£ = 0, . . . , k — 1, m = 0, . . . , p — 1) correspond to the standard coefficients of the k-step 
Adams-Bashforth scheme (IS) and the coefficients of the LTS-ABk(p) scheme defined in (21), respectively. 



Proof. The identity in (29) was verified by computer algebra for all k < 20 and p < 1000, which is sufficient 
for all practical purposes. It probably holds for all values of k and p. □ 

We are now ready to establish the accuracy of the LTS-ABfc(p) scheme (Algorithm 3.1). 

Theorem 1. The local time-stepping method LTS-ABk(p) is consistent of order k. 

Proof. For simplicity, we restrict ourselves to the cases k = 2, 3 and 4, as the extension to the general case 
k > 4 is straightforward but cumbersome. 

To prove the second-order consistency of LTS-AB2(p), we need to show that the local error y n +i — y(i n +i) 
is 0(h 3 ). Since (26) with k — 2 and At — At/p holds for all m = 0, . . . ,p — 1, we have that 



At, 



o-l 



p-1 



y«+i = yi = yo + —^i 1 ~ p ) ( ^2 ^™^y™ + ^ Pm,iyn-i 



At 

p 



\m— m= 

p-2 



BP a y(p-i)/ P + ("o + ai)y m / P + "iy-1/p 



Next, we use that y :=y n as well as (28) and (29) with k — 2. This yields 

/ 



y„+i = Yn + AtB(I - P) (a y n + aiy„_i) + —BP 

P 



a oy(p-i)/ P + ^ y-m/p + "iy-1/p 

V m=0 J 



(30) 



For r = and k — 2, we find from (21) and 7j(0) = 0, j > 1, that 

y'(0) = B(I - P)y„ + BPy(0) = By„ = y'(t n ) . 
Thus, we may expand y m / p , m = — 1, . . . ,p — 1, in Taylor series as 

Vm/p = y(0) + - Aty'(O) + 0(At 2 ) = y„ + - Aty'(t n ) + 0(At 2 ) . 



P 



P 



(31) 



We now insert (31) into (30), replace y n ~i by its Taylor expansion, use (28) with k = 2 and obtain after 
some simplifications 

y„+i = y„ + AtB(I - P) (y„ - ai Aty'(t n ) + <D{At 2 )) + AtBP (y„ +C X Aty'(i„) + G(At 2 )) , 
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where 



(p-l)(p-2) 
cw - !) H n "i 



With the coefficients of the classical two-step Adams-Bashforth scheme ao — 3/2 and a\ = —1/2, we note 
that 

ai+Ci = i [-p 2 +3(p-l) + (p-l)(p-2) + l] -0. 



By differentiating (15), we thus find 



At 2 



y n +i =y n + At B (y„ - ai Aty'(t n ) + 0(At 2 )) = y„ + Aty'(i n ) + — y"(t„) + 0(At 3 
which yields a local error of order 0(At 3 ). 



For fc = 3, we need to prove a local error of order 0(At 4 ). Similarly to the derivation of ( 30 ) , we now 



find 



y n +i = yn + At B(I - P) (a y„ + aiy„_i + a 2 y„_ 2 ) 
Ai, 



P 



/ P-3 

BP I a y( p _i)/p + (ao + «i)y(p-2)/p + ^2 

\ m=0 



Ym/p + (ai + a 2 )y_i/ p + a 2 y_ 2/p 



(32) 



Again, from (21) we have y'(0) = y'(t„). By differentiation of (21) and using Taylor expansions we also find 
y"(0) = B(I - P) XI %(0)V j Yn + BPy'(O) = B(I - P)i- f|y„ - 2y„_! + Jy„_ 2 ) + BPy'(t n ) 



3=0 
1 



B(I - P)— (Aty'(t„) + 0(At 2 )) + BPy'(f n ) = y"(t„) + O(At) . 



At 

Thus, we deduce that 



Aty'(t r> 



p ' "' p 2 2 

By following similar arguments as for k = 2, we now obtain 



y"(t„) + 0(At 3 ) , m = -2 J ...,p-l. 



At 2 



y„+i = y„ + At B(I - P) ( y n - (ax + 2a 2 ) Aty'{t n ) + (ai + 4a 2 ) — y"(t n ) + 0(At 3 ) 

At 2 

At BP ( y„ + d Aty'(t„) + C 2 — y"(t„) + 0(At 3 ; 



where 



Ci := 



Co := 



a (p - 1) + (a + a x )(p - 2) + ^ — ^ — ^ - (ai + a 2 ) - 2a 2 



a (p ~ I) 2 + (ao + ai)(p - 2) 2 + 



(p-2)(p-3)(2p-5) 
6 



+ (ai + a 2 ) + 4a 2 



(33) 



With the coefficients of the classical three-step Adams-Bashforth scheme ao = 23/12, ol\ = —16/12 and 
a 2 = 5/12, we can easily verify that 

a 1 + 2a 2 + C\ = , —a\ — 4a 2 + C 2 = , — ai — 2a 2 = \- , ai + 4a 2 = - . 



Hence ( 33 ) reduces to 



At 2 At 3 
y n+1 = y n + Aty'(t„) + — y"(t„) + — y"'(t„) + C(At 4 ) . 
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which completes the proof for the LTS-AB3(p) scheme. 

The case k = 4 follows from similar arguments. For k = 4, we find from (21 1 with r = that y'(0) 
y'(t n ) and by using Taylor expansions that 



y"(0) = B(I-P) 



1 /ll 



At V 6 
1 



3y n -i 



5 1 

;yn-2 — TiYnS 



■BPy'(O)=y"(i n )+0(At), 



y"'(0) = B(I - P)^ - 5 y«-i + 4 y«-2 - Yn-a) + BPy"(0) = y"'{t n ) + 0{At) 



This yields 



Ym/p =Yn-\ Aty (t n ) 

' p 



p 2 2 



y"(t n ) 



TO 3 At 3 

p 3 6 



y"'(t n )+0(At 4 ), m=-3,...,p-l. (34) 



Next, we insert (34 1 into (26 1 with fc = 4 and At = At/p, replace y n _i, y n -2 and y n -3 by their respective 
Taylor expansions, use (28 1 and (291 with k = 4, and find after further simplifications 



At 2 At 3 At 4 

y n+1 = y n + Aty'(t n ) + — y"(t n ) + — y"'(t n ) + — y""(t n ) + 0(At 5 ) . 



which completes the proof for the LTS-AB4(p) scheme. 



□ 



4. Numerical experiments 

Here we present numerical experiments that illustrate the stability properties of the above LTS methods, 
validate their expected order of convergence and demonstrate their usefulness in the presence of complex 
geometry. First, we consider a simple one-dimensional test problem to study the stability of the different 
LTS schemes presented above and to show that they yield the expected overall rate of convergence when 
combined with a spatial finite clement discretization of comparable accuracy, independently of the number 
of local time-steps p used in the fine region. Then, we illustrate the versatility of our LTS schemes by 
simulating the propagation of a plane wave in a rectangular domain adjacent to a roof mounted antenna. 

4-1- Stability 

We consider the one-dimensional homogeneous damped wave equation (JlJ with constant wave speed 
c = 1 and damping coefficient a — 0.1 on the interval f2 = [0 , 6]. Next, we divide £1 into three equal parts. 
The left and right intervals, [0, 2] and [4, 6], respectively, are discretized with an equidistant mesh of size 
^coarse ^ wnereas interval flf — [2,4] is discretized with an equidistant mesh of size /i fine = h com ' sc /p - 
see Fig. [2j Hence, the two outer intervals correspond to the coarse region whereas the inner interval [2 , 4] 
to the refined region. 

For every time-step At, we shall take p steps of size At = At/p inside the refined region £lf. In the 
absence of local refinement, i.e., p — 1, the mesh is equidistant throughout ft. Then, the LTS-ABfc(p) 
algorithm reduces to the standard k-step Adams-Bashforth method and we denote by AtABk the largest 
time-step allowed. For p > 2, we let At p denote the maximal time-step permitted in the LTS-ABfc(p) 
algorithm; clearly, we always have At p < AtABk- When At p = AtABk, the LTS algorithm imposes no 
further restriction on At; we then call the CFL condition of the LTS scheme optimal. We shall now evaluate 
numerically the CFL condition of the various LTS schemes from Section 3. To do so, we proceed as follows: 

1. Set AtABk to the maximal At allowed for the equidistant mesh of mesh size h coalsc ; 

2. refine the equidistant mesh p times inside f2/; 

3. determine the maximal time-step At p allowed by the LTS-ABfc(p) method on the locally refined mesh 
and compare At p with AtABk- 
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iiiiiiiiii 



i iiiiiiiiiiiiiiiiiiii 



Figure 2: One-dimensional example: the computational domain Q = [0,6] with the refined region Qf = [2,4]. 



First, we consider a? 1 continuous FE discretization with mass lumping and combine it with the second- 
order Adams-Bashforth scheme. We choose h com ' sc = 0.1 and rewrite the two-step Adams-Bashforth method 
as the one-step method 

vii\ / v \ / I + ilAfB — ^-B \ 



l = C AB2 I I , WB2 . 

It is stable if p{Cab2) < lj where p{Cab2) denotes the spectral radius of the matrix Cab2 [32] • Progressively 
increasing At while monitoring p(Cab2), we find that the maximal time-step allowed is AtAB2 — 0.0106 for 
p = 1. Next, we refine by a factor p = 2 those elements that lie inside the interval [2, 4], that is h — 0.05, 
and set to one all corresponding entries in P. Hence, for every time-step At, we shall take two steps of size 
At = At/2 inside the refined region with the second-order time-stepping scheme LTS-AB2(2). To determine 
the stability of the LTS-AB2(2) scheme, we also rewrite it as a one-step method: 

ClTS-AB2(2) Y-l/2 j C-LTS-AB2(2) 





Cu 


C12 


C13 




C22 


C23 


I 









with 



C^ I+ ^B-I it BP^(f) 3 ,BP )2 + f(^) 2 BPB, 

C " = -J AfBP -i(f) 2(BP)2 ' 

i- + l*HP + 5(-) 2 ( BP,-|(^BPB, 

C 2 i = I+ -AB+-At BP , C 22 = --AfBP, C 23 = -- AB+ - AtBP . 

8 8 4 6 8 8 

The LTS-AB2(2) scheme is stable if p(C LTS _ AB2 (2)) < 1) where p{Clts-ab2(2)) denotes the spectral radius 
of the matrix C LTS _ AB2{2 ). 

To determine the range of values At for which the LTS-AB2(2) scheme is stable, we compute the spectral 
radius of Clts-ab2(2) f° r varying At. As shown in the right frame of Fig. [3l the spectral radius transgresses 
the stability threshold at one for a time-step about 80% of At A B2- Thus, for all time-steps At < 0.8 AtAB2, 
the LTS-AB2(2) scheme is stable; yet its CFL condition is not optimal. 

Next, we consider the third-order Adams-Bashforth scheme and combine it with a V 2 continuous FE 
discretization with mass lumping. We choose /j coarso = 0.2, which yields the maximal time-step AtAB3 = 
0.029 for p = 1. Again, we refine by a factor p — 2 those elements that lie inside the interval [2,4], that is 
^fmc _ q ^ Hence, for every time-step At, we shall take two steps of size At = At/2 in the refined region 
with the third-order LTS-AB3(2) scheme. 

To study its stability properties, we again reformulate both the classical AB3 and the LTS-AB3(2) 
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Figure 3: The classical AB2 and the LTS-AB2(2) schemes combined with V 1 continuous FE: the spectral radius of Cab2 (left) 
and Cj/rS— j4.B2(2) (fight) is shown for varying At/ AtAB2- 



0.5 0.6 0.7 0.0 0.9 

At/At 



Figure 4: The classical AB3 and the LTS-AB3(2) schemes combined with "P 2 continuous FE: the spectral radius of C^g^ (left) 
and C iTS _ AS3 (2) (right) is shown for varying At/AtAB3- 



schemes as one-step methods: 




J AB3 




( y«+i \ 
yi/2 

\ Yn-1 / 



J LTS-AB3{2) 



( y n 

y-1/2 
y«-i 
\ y n -2 



and compute the spectral radius of Cab3 and Clts-AB3(2) f° r varying At. In the right frame of Fig. |4j 
we observe that the spectral radius of Clts-ab3(2) lies below one for all time-steps At < AtAB3- Hence, 
the LTS-AB3(2) scheme is stable up to the maximal time-step allowed by the standard third-order Adams- 
Bashforth method on an equidistant mesh; therefore its CFL condition is optimal. 

Finally, to study the stability of the fourth-order LTS scheme, we consider a V z continuous FE dis- 
cretization with mass lumping. Again, we choose /j coarsc = 0.2, which now yields the maximal time-step 
At aba = 0.0099 for p = 1 and refine by a factor p = 2 those elements that lie inside the interval [2, 4]. Hence, 
for every time-step At, we shall use the fourth-order time-stepping scheme LTS-AB4(2) with At = At/2 in 
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1.03 



Figure 5: The classical AB4 and the LTS-AB4(2) schemes combined with "P 3 continuous FE: the spectral radius of Gab4 (left) 
and C LT g_A Bi (2) ( r ight) is shown for varying At/AtABi- 



Figure 6: The fc-th order LTS-ABfc(2) scheme combined with IP-DG V ^elements: the spectral radius of CLTS—ABk(2) with 
k = 2 (left) and k = 3 (right) is shown for varying At/ AtABk- 



the refined region. After refomulating the AB4 and the LTS-AB4(2) schemes as one-step methods 



/ y«+i \ 

y«-i 
V y«-2 J 



= c 



AB4 



( y» \ 

yn-l 
yn-2 
\ yn-3 J 



( y n +i \ 
yi/2 
y n 

y-1/2 
y n -i 

V y-3/2 J 



= c 



LTS~AB4(2) 



( y w 

y-1/2 
y n -i 
y-3/2 
y«-2 
V y«-3 



we compute the spectral radius of Caba and Clts-ab4(2) for varying At/ AtABi- As shown in the right 
frame of Fig. [5j the spectral radius of Clts-ab4{2) lies below one for all time-steps At < AtAB4- Thus, the 
LTS-AB4(2) scheme is stable up to the optimal time-step. 

So far, we have restricted the detailed numerical stability study of the LTS schemes to standard con- 
tinuous FE discretizations, where the accuracy of the spatial discretization matches that of the time dis- 
cretization. We obtain similar stability results, when the LTS-ABfc schemes are combined with a fc-th order 
IP-DG or nodal DG discretization in space. For instance, as shown in the left frame of Fig. |6j the largest 
time-step allowed by the LTS-AB2(2) scheme when combined with IP-DG 'P 1 -elements again is about 80% 
of AtAB2- Similarly, the spectral radius of C LT s_ab3(2)> shown in the right frame of Fig. [6] for varying 
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At/ AtAB3, reveals that the LTS-AB3(2) scheme, when combined with IP-DG T^-elements, is again stable 
for the optimal time-step AtAB3- 

Remark 1. In summary, our numerical experiments show for a spatial discretization with standard contin- 
uous, IP-DG or nodal DG finite elements: 

• the maximal time-step At p allowed by the LTS-AB2(p) scheme is about 80 % of the optimal time-step 
Atj±B2 for all values of h, p and a; 

• the CFL stability condition of the LTS-AB3(p) and LTS-ABA(p) schemes is optimal for all h, p and 
a. 



4-2. Convergence 

We consider the one-dimensional homogeneous model problems (JlJ and Q with constant wave speed 
c = 1 and damping coefficient a — 0.1 on the interval Q — [0 , 6]. The initial conditions are chosen to yield 
the exact solution 

u(x, t) = = sin(7ra:) sin ( - \J Att 2 — a 2 

^ ' VAir 2 - a 2 V ' V 2 

On 

v(x,t) = — (x,t) , vr(x,t) = -Vu(x,t) . 

Again, we divide f2 into three equal parts. The left and right intervals, [0, 2] and [4, 6], respectively, are dis- 
cretized with an equidistant mesh of size h c ° msc , whereas the interval [2, 4] is discretized with an equidistant 
mesh of size h Rnc = ft, coarsc /p. Hence, the two outer intervals correspond to the coarse region and the inner 
interval [2 , 4] to the refined region. The first k — 1 time-steps of each LTS-ABfc(p) scheme are initialized by 
using the exact solution. 




(a) continuous FE (h = 0.02, 0.01, 0.005, 0.0025) (b) IP-DG (h = 0.02, 0.01, 0.005, 0.0025) 



— P- 2 
^p-5 

— P-7 




(c) nodal DG (h = 0.02, 0.01, 0.005, 0.0025) 
Figure 7: LTS-AB2(p) error vs. h = h coaISC for V 1 finite elements with p = 2, 5, 7. 
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First, we consider a V 1 continuous FE discretization with mass lumping and a sequence of increasingly 
finer meshes. For every time-step At, we shall take p > 2 local steps of size At = At/p in the refined region, 
with the second-order time-stepping scheme LTS-AB2(p). According to our previous results on stability 
from Section 4.1, we set At = 0.8 • AtAB2\ note that AtAB2 depends on the mesh size. As we systematically 
reduce the global mesh size /i coarsc ; while simultaneously reducing At, we monitor the L 2 space-time error 
in the numerical solution \\u(-,T) — u h (-,T)\\ L 2^ at the final time T = 10. In frame (a) of Fig. [TJ the 
numerical error is shown vs. the mesh size h — h coa - lsc : regardless of the number of local time-steps p = 2, 5 
or 7, the numerical method converges with order two. 

We now repeat the same experiment with the IP-DG (a = 5 in ([6])) and the nodal DG discretizations 
with T-^-elements. As shown in frames (b) and (c) of Fig. [7J the LTS-AB2(p) method again yields overall 
second-order convergence independently of p. 

Next, we consider the third-order LTS-AB3(p) scheme and combine it with any one of the three FE 
discretizations with 7-' 2 -elements. Thus, we expect all numerical schemes to exhibit overall third-order 
convergence with respect to the L 2 -norm. We set At = AtAB3, the largest possible time-step allowed by the 
third-order Adams-Bashforth approach on an equidistant mesh with h = h coarsc . In Fig. [8] we display the 
space-time L 2 -errors of the numerical solutions at T = 10 for a sequence of meshes and different values of p. 
The continuous FEM with mass lumping, the IP-DG method (with a = 12) and the nodal DG discretization 
all yield the expected third-order convergence. 

Finally, to demonstrate the order of convergence of the fourth-order LTS-AB4(p) scheme, we consider 
again the continuous FE or the two DG discretizations with T-^-elements. Here, we set the penalty parameter 
a — 20 for the IP-DG method and let At — At aba, the corresponding largest possible time-step allowed by 
the Adams-Bashforth approach of order four on an equidistant mesh with h = /j coarse . We monitor the I? 
space-time error in the numerical solution at T = 10 for the sequence of meshes and different values of p. 
Again, the numerical results shown in Fig. [9] corroborate the expected fourth-order convergence. 
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(a) continuous FE (h = 0.2, 0.1, 0.05, 0.025) (b) IP-DG (h = 0.2, 0.1, 0.05, 0.025) 




(c) nodal DG (h = 0.02, 0.01, 0.005, 0.0025) 
Figure 9: LTS-AB4(p) error vs. h = h coarsc for V 3 finite elements with p = 2, 5, 7. 



Remark 2. We obtain similar convergence results for other values of p and a. In summary, we observe 
the expected convergence of order k for the LTS-ABk(p) schemes, regardless of the spatial FE discretization 
and independently of the number of local time-steps p and the damping coefficient a. 

4-3. Two-dimensional example 

To illustrate the usefulness of the LTS approach, we consider in Q, a rectangular domain of size [0, 3] x 
[0, 1] adjacent to a roof mounted antenna of thickness 0.01, as shown in Fig. 10 We set the constant wave 
speed c — 1 and the constant damping coefficient a = 0.1. On the boundary of f2 we impose homogeneous 
Neumann conditions and choose as initial conditions the vertical Gaussian plane wave 

Uq(x, y) = exp (-(x - x f/S 2 ) , for all (x, y) £ D, , 
v (x, y) = , for all (x, y) £ Q, , 

centered about xq = 1.3 and of width 5 = 0.01. 

For the spatial discretization we opt for V 2 continuous finite elements with mass lumping. First, fl is 
discretized with triangles of maximal size /j coarso = 0.05. However, such coarse triangles do not resolve 
the small geometric features of the antenna, which require h w /i coarsG /7, as shown in Fig. 10 Then, we 



h coarsc /7, as shown in Fig. 

successively refine the entire mesh three times, each time splitting every triangle into four. Since the initial 
mesh in Q is unstructured, the boundary between the fine and the coarse mesh is not well-defined. Here 
we let the fine mesh correspond to all triangles with h < 0.6h coaisc in size, i.e. the darker (green) triangles 



in Fig. 10 The corresponding degrees of freedom in the finite element solution are then selected merely by 
setting to one the corresponding diagonal entries of the matrix P (see Section 3.2). 

For the time discretization, we choose the third-order LTS-AB3(7) time-stepping scheme with p = 7, 
which for every time-step At takes seven local time-steps inside the refined region that surrounds the antenna 



(see Fig. 10 1. Thus, the numerical method is third-order accurate in both space and time under the CFL 
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Figure 10: The initial triangular mesh (left); zoom on the "fine" mesh indicated by the darker (green) triangles (right). 



Figure 11: Two-dimensional example: the solution is shown at times t =0.4, 0.55, 0.71, 0.82, 1 and 1.2. 



condition At = 0.07 /i coalsc , determined experimentally. If instead the same (global) time-step At was used 
everywhere inside tt, it would have to be about seven times smaller than necessary in most of tt. As a 
starting procedure, we employ a standard fourth-order Runge-Kutta scheme. 



In Fig. 11 snapshots of the numerical solution are shown at different times. The vertical Gaussian pulse 
initiates two plane waves propagating in opposite directions. At t — 0.5, the right-moving wave impinges 
first on the antenna and than on the tip of roof. Multiple reflections occur as the lower part of the wave 
bounces back. Meanwhile, the upper part of the plane wave has proceeded to the right without any spurious 
reflection between the coarse and the refined regions. 



5. Conclusion 

Starting from classical Adams-Bashforth methods, we have presented explicit local time-stepping (LTS) 
schemes for damped wave equations, which permit arbitrarily small time-steps precisely where the smallest 
elements in the mesh are located. Thus, when combined with a spatial finite element discretization with 
an essentially diagonal mass matrix, the resulting time-marching schemes are fully explicit. The general 
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fc-th order LTS scheme, denoted by LTS-ABfc, is given by the LTS-ABfc(p) Algorithm in Section 3.2. As 
shown in Section 3.4, it is fc-th order accurate in time. Thus when combined with a spatial finite element 
discretization of order fc— 1, the resulting numerical scheme yields optimal fc-th order space-time convergence 
in the L 2 -norm. The derivation of the LTS-ABfc (p) algorithm immediately generalizes to the inhomogeneous 
case with nonzero external forcing. 

The LTS-ABfc(p) scheme has been combined with three distinct finite element discretizations: standard 
H 1 -conforming finite elements, an IP-DG formulation, and nodal DG finite elements. In all cases, our 
numerical results demonstrate that the resulting LTS-ABfc schemes of order k > 3 have optimal CFL 
stability properties regardless of the mesh size h, the global to local step-size ratio p, or the dissipation 
a. For k = 2, however, the CFL condition is sub-optimal as the maximal time-step allowed by the LTS- 
AB2 scheme is only about 80% of that permitted by the standard second-order Adams-Bashforth scheme 
on an equidistant mesh. In contrast to the energy conserving LTS methods that were developed for the 
(undamped) wave equation in [21] . the LTS-ABfc methods with k > 3 always achieve the optimal CFL 
condition without overlap of the fine and the coarse region. 

Since the LTS methods presented here are truly explicit, their parallel implementation is straightforward. 
Let At denote the time-step imposed by the CFL condition in the coarser part of the mesh. Then, during 
every (global) time-step At, each local time-step of size At/p inside the fine region of the mesh, simply 
corresponds to sparse matrix- vector multiplications that only involve the degrees of freedom associated with 
the fine region of the mesh. Those "fine" degrees of freedom can be selected individually and without any 
restriction by setting the corresponding entries in the diagonal projection matrix P to one; in particular, no 
adjacency or coherence in the numbering of the degrees of freedom is assumed. Hence the implementation 
is straightforward and requires no special data structures. 

In the presence of multi-level mesh refinement, each local time-step in the fine region can itself include 
further local time-steps inside a smaller subregion with an even higher degree of local mesh refinement. The 
explicit local time-stepping schemes developed here for the scalar damped wave equation immediately apply 
to other damped wave equations, such as in electromagnetics or elasticity; in fact, they can be used for 
general linear first-order hyperbolic systems. 
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